#include <cmath>
#include <limits>
#include <boost\numeric\interval.hpp>
#include <gmpxx.h>
#include "predicates.h"

int left_turn(const point& a, const point& b, const point& c)
{
    double left_turn = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);

    if (eps * (std::abs((b.x - a.x) * (c.y - a.y)) + std::abs((b.y - a.y) * (c.x - a.x))) < left_turn)
    {
        return 1;
    }

    if (-(eps * (std::abs((b.x - a.x) * (c.y - a.y)) + std::abs((b.y - a.y) * (c.x - a.x)))) > left_turn)
    {
        return -1;
    }

    return exact_left_turn(a, b, c);
}

void value_left_turn(const point& a, const point& b, const point& c, mpq_t& result)
{
    mpq_t bx;
    mpq_init(bx);
    mpq_set_d(bx, b.x);

    mpq_t ax;
    mpq_init(ax);
    mpq_set_d(ax, a.x);

    mpq_t cy;
    mpq_init(cy);
    mpq_set_d(cy, c.y);

    mpq_t ay;
    mpq_init(ay);
    mpq_set_d(ay, a.y);

    mpq_t by;
    mpq_init(by);
    mpq_set_d(by, b.y);

    mpq_t cx;
    mpq_init(cx);
    mpq_set_d(cx, c.x);

    mpq_t tmp1;
    mpq_init(tmp1);
    mpq_sub(tmp1, bx, ax);

    mpq_t tmp2;
    mpq_init(tmp2);
    mpq_sub(tmp2, cy, ay);

    mpq_mul(result, tmp1, tmp2);

    mpq_sub(tmp1, by, ay);
    mpq_sub(tmp2, cx, ax);
    mpq_mul(tmp1, tmp1, tmp2);
    mpq_sub(result, result, tmp1);

    mpq_clear(ax); mpq_clear(bx); mpq_clear(cx); mpq_clear(ay); mpq_clear(by); mpq_clear(cy); mpq_clear(tmp1); mpq_clear(tmp2);
}

int exact_left_turn(const point& a, const point& b, const point& c)
{
    mpq_t value;
    mpq_init(value);
    value_left_turn(a, b, c, value);

    int sign = mpq_sgn(value); mpq_clear(value);

    return sign;
}

inline double sqrradius(const point& a)
{
    return sqrhypot(a, point(0, 0));
}

template <class T>
void det_multiplied_by_second_det(const point& p, const circle& circ, T& result)
{
    point a = circ.boundaries[0]; point b = circ.boundaries[1]; point c = circ.boundaries[2];

    T r_a(sqrradius(a));
    T r_b(sqrradius(b));
    T r_c(sqrradius(c));

    T a_x(a.x), b_x(b.x), c_x(c.x), a_y(a.y), b_y(b.y), c_y(c.y);

    T first_det = (b_x - a_x) * (c_y - a_y) - (b_y - a_y) * (c_x - a_x);
    T second_det = (r_b - r_a) * (c_y - a_y) - (b_y - a_y) * (r_c - r_a);
    T third_det = (r_b - r_a) * (c_x - a_x) - (b_x - a_x) * (r_c - r_a);
    T fourth_det = r_a * (b_x * c_y - c_x * b_y) + r_b * (c_x * a_y - a_x * c_y) + r_c * (a_x * b_y - b_x * a_y);

    T r(sqrradius(p));
    T x(p.x);
    T y(p.y);

    T det = r * first_det - x * second_det + y * third_det - fourth_det;
    
    result = det * first_det;
}

bool point_in_circle(const point& p, const circle& circ)
{
    switch(circ.boundaries.size())
    {
    case 1:
        return p == circ.boundaries[0];
    case 2:
    {
        double r_bp = sqrhypot(circ.boundaries[1], p);
        double r_ap = sqrhypot(circ.boundaries[0], p);
        double r_ab = sqrhypot(circ.boundaries[0], circ.boundaries[1]);

        return r_ap + r_bp <= r_ab;
    }
    case 3:
    {
        boost::numeric::interval<double> result;
        det_multiplied_by_second_det(p, circ, result);

        if ((result.lower() <= 0) && (result.upper() >= 0))
        {
            mpq_class res;
            det_multiplied_by_second_det(p, circ, res);
            
            return res <= 0;
        } else
        {
            return result.upper() <= 0;
        }
    }
    default:
        return false;
    }
}

int further_from_line_than_old_uppermost(const point& a, const point& b, const point& c, const point& d)
{
    double result = (b.x - a.x) * (c.y - d.y) - (b.y - a.y) * (c.x - d.x);

    if (eps * (std::abs((b.x - a.x) * (c.y - d.y)) + std::abs((b.y - a.y) * (c.x - d.x))) < result)
    {
        return 1;
    }

    if (-(eps * (std::abs((b.x - a.x) * (c.y - d.y)) + std::abs((b.y - a.y) * (c.x - d.x)))) > result)
    {
        return -1;
    }
    
    return exact_further_from_line_than_old_uppermost(a, b, c, d);
}


int exact_further_from_line_than_old_uppermost(const point& a, const point& b, const point& c, const point& d)
{
    mpq_t bx;
    mpq_init(bx);
    mpq_set_d(bx, b.x);

    mpq_t ax;
    mpq_init(ax);
    mpq_set_d(ax, a.x);

    mpq_t cy;
    mpq_init(cy);
    mpq_set_d(cy, c.y);
    
    mpq_t dy;
    mpq_init(dy);
    mpq_set_d(dy, d.y);

    mpq_t by;
    mpq_init(by);
    mpq_set_d(by, b.y);

    mpq_t ay;
    mpq_init(ay);
    mpq_set_d(ay, a.y);

    mpq_t cx;
    mpq_init(cx);
    mpq_set_d(cx, c.x);

    mpq_t dx;
    mpq_init(dx);
    mpq_set_d(dx, d.x);

    mpq_t tmp1;
    mpq_init(tmp1);
    mpq_sub(tmp1, bx, ax);
    //tmp1 = bx - ax

    mpq_t tmp2;
    mpq_init(tmp2);
    mpq_sub(tmp2, cy, dy);
    //tmp2 = cy - dy

    mpq_t result;
    mpq_init(result);
    mpq_mul(result, tmp1, tmp2);
    //result = (bx - ax) * (cy - dy)

    mpq_sub(tmp1, by, ay);
    //tmp1 = by - ay
    mpq_sub(tmp2, cx, dx);
    //tmp1 = cx - dx
    mpq_mul(tmp1, tmp1, tmp2);
    //tmp1 = (by - ay) * (cx - dx)
    mpq_sub(result, result, tmp1);
    //result = (bx - ax) * (cy - dy) - (by - ay) * (cx - dx)

    mpq_clear(ax); mpq_clear(bx); mpq_clear(cx);  mpq_clear(dx); mpq_clear(ay); mpq_clear(by); mpq_clear(cy); mpq_clear(dy); mpq_clear(tmp1); mpq_clear(tmp2);

    int sign = mpq_sgn(result); mpq_clear(result);

    return mpq_sgn(result);
}